end
solve_ICCG （ 7/7 ） : METHOD=1
!C
!C +---+
!C | ALPHA= RHO / {p}{q} |
!C +---+
!C===
C1= 0.d0 do i= 1, N
C1= C1 + W(i,P)*W(i,Q) enddo
ALPHA= RHO / C1
!C===
!C
!C +---+
!C | {x}= {x} + ALPHA*{p} |
!C | {r}= {r} - ALPHA*{q} |
!C +---+
!C===
do i= 1, N
X(i) = X(i) + ALPHA * W(i,P) W(i,R)= W(i,R) - ALPHA * W(i,Q) enddo
DNRM2= 0.d0 do i= 1, N
DNRM2= DNRM2 + W(i,R)**2 enddo
!C===
ERR = dsqrt(DNRM2/BNRM2) if (ERR .lt. EPS) then
IER = 0 goto 900 else
RHO1 = RHO endif
enddo IER = 1
Compute r
(0)= b-[A]x
(0)for i= 1, 2, …
solve [M]z
(i-1)= r
(i-1)
i-1= r
(i-1)z
(i-1)if i=1
p
(1)= z
(0)else
i-1=
i-1/
i-2p
(i)= z
(i-1)+
i-1p
(i)endif
q
(i)= [A]p
(i)
i=
i-1/p
(i)q
(i)x
(i)= x
(i-1)+
ip
(i)r
(i)= r
(i-1)-
iq
(i)check convergence |r|
r= b-[A]x end
DNRM2=|r|
2BNRM2=|b|
2ERR= |r|/|b|
solve_ICCG2 （ 1/3 ） : METHOD=2
!C
!C***
!C*** module solver_ICCG2
!C***
!
module solver_ICCG2 contains
!C
!C*** solve_ICCG2
!C
subroutine solve_ICCG2 &
& ( N, NPL, NPU, indexL, itemL, indexU, itemU, D, B, X, &
& AL, AU, EPS, ITR, IER) implicit REAL*8 (A-H,O-Z)
real(kind=8), dimension(N) :: D real(kind=8), dimension(N) :: B real(kind=8), dimension(N) :: X real(kind=8), dimension(NPL) :: AL real(kind=8), dimension(NPU) :: AU
integer, dimension(0:N) :: indexL, indexU integer, dimension(NPL):: itemL
integer, dimension(NPU):: itemU
real(kind=8), dimension(:,:), allocatable :: W integer, parameter :: R= 1
integer, parameter :: Z= 2 integer, parameter :: Q= 2 integer, parameter :: P= 3 integer, parameter :: DD= 4
real(kind=8), dimension(:), allocatable :: ALlu0, AUlu0 real(kind=8), dimension(:), allocatable :: Dlu0
Compute r
(0)= b-[A]x
(0)for i= 1, 2, …
solve [M]z
(i-1)= r
(i-1)
i-1= r
(i-1)z
(i-1)if i=1
p
(1)= z
(0)else
i-1=
i-1/
i-2p
(i)= z
(i-1)+
i-1p
(i-1)endif
q
(i)= [A]p
(i)
i=
i-1/p
(i)q
(i)x
(i)= x
(i-1)+
ip
(i)r
(i)= r
(i-1)-
iq
(i)check convergence |r|
end
solve_ICCG2 （ 2/3 ） : METHOD=2
!C
!C +---+
!C | INIT |
!C +---+
!C===
allocate (W(N,4)) do i= 1, N
X(i) = 0.d0 W(i,1)= 0.0D0 W(i,2)= 0.0D0 W(i,3)= 0.0D0 W(i,4)= 0.0D0 enddo
call FORM_ILU0
!C===
Dlu0, ALlu0
，AUlu0
にはILU(0)
分解された対角，下三角，上三角成分が入る（行列
[M]
）。不完全修正コレスキー分解：正確には不完全修正 LU 分解 solver_ICCG2.f に附属
contains
!C
!C***
!C*** FORM_ILU0
!C***
!C
!C form ILU(0) matrix
!C
subroutine FORM_ILU0 implicit REAL*8 (A-H,O-Z)
integer, dimension(:), allocatable :: IW1 , IW2 integer, dimension(:), allocatable :: IWsL, IWsU real (kind=8):: RHS_Aij, DkINV, Aik, Akj
!C
!C +---+
!C | INIT. |
!C +---+
!C===
allocate (ALlu0(NPL), AUlu0(NPU)) allocate (Dlu0(N))
do i= 1, N Dlu0(i)= D(i) do k= 1, INL(i)
ALlu0(k,i)= AL(k,i) enddo
do k= 1, INU(i)
AUlu0(k,i)= AU(k,i) enddo
enddo
!C===
「
Dlu0,ALlu0,AUlu0
」初期値として，「D,AL,AU」の値を代入する。
jk k j
k ik ij
ij
a l d l
l
1
1
n i 1 , 2 , ,
1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
Dlu0, ALlu0，AUlu0にはILU(0)分解
された対角，下三角，上三角成分が 入る（行列[M]
）。!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
do i= 1, N do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo endif enddo enddo
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
do i= 1, N
do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo
endif
enddo
enddo
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
do i= 1, N
do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo
endif
enddo
enddo
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
do i= 1, N
do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo endif enddo enddo
j=i
jk k j
k ik ij
ij
a l d l
l
1
1
n i
1,2,,1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
jk k j
k ik ij
ij
a l d l
l
1
1
n i
1,2,,1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
do i= 1, N
do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo endif enddo enddo
j<i
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
do i= 1, N
do k= 1, i-1
if (A(i,k) is non-zero) then do j= k+1, N
if (A(i,j) is non-zero) then
A(i,j)= A(i,j) &
-A(i,k)*(A(k,k))
-1*A(k,j) endif
enddo endif enddo enddo
j>i
jk k j
k ik ij
ij
a l d l
l
1
1
n i
1,2,,1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
実はこのIF文の中は通らない
（理由は後述）したがって，
ALlu0= AL AUlu0= AU
!C
!C +---+
!C | ILU(0) factorization |
!C +---+
!C===
allocate (IW1(N) , IW2(N)) IW1=0
IW2= 0 do i= 1, N
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= k0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= k0
enddo
do icon= indexL(i-1)+1, indexL(i) k= itemL(icon)
D11= Dlu0(k) DkINV= 1.d0/D11
Aik= ALlu0(icon)
do kcon= indexU(k-1)+1, indexU(k) j= itemU(kcon)
if (j.eq.i) then
Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj Dlu0(i)= Dlu0(i) - RHS_Aij endif
if (j.lt.i .and. IW1(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW1(j)
ALlu0(ij0)= ALlu0(ij0) - RHS_Aij endif
if (j.gt.i .and. IW2(j).ne.0) then Akj= AUlu0(kcon)
RHS_Aij= Aik * DkINV * Akj ij0 = IW2(j)
AUlu0(ij0)= AUlu0(ij0) - RHS_Aij endif
enddo enddo
do k0= indexL(i-1)+1, indexL(i) IW1(itemL(k0))= 0
enddo
do k0= indexU(i-1)+1, indexU(i) IW2(itemU(k0))= 0
enddo enddo do i= 1, N
Dlu0(i)= 1.d0 / Dlu0(i) enddo
deallocate (IW1, IW2)
!C===
end subroutine FORM_ILU0
j>i
jk k j
k ik ij
ij
a l d l
l
1
1
n i
1,2,,1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
j<i
i
k
ある要素「i
（○）」に接続する下三角成 分「k
（■○）」の上三角成分「j
（■）」が：（
1
）「
i
」自身である場合，Dlu
（■）更新j=i
j
jk k j
k ik ij
ij
a l d l
l
1
1
n i 1 , 2 , ,
1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
i
k
ある要素「i
（○）」に接続する下三角成 分「k
（■○）」の上三角成分「j
（■）」が：（
2
）「
i
」の下三角成分である場合ALlu0(i-j)
（■）更新j<i j
jk k j
k ik ij
ij
a l d l
l
1
1
n i 1 , 2 , ,
1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
i
k
ある要素「i
（○）」に接続する下三角成 分「k
（■○）」の上三角成分「j
（■）」が：（
3
）「
i
」の上三角成分である場合AUlu0(i-j)
（■）更新実際は（
2
），（3
）に該当する場合は無しj>i
j
jk k j
k ik ij
ij
a l d l
l
1
1
n i 1 , 2 , ,
1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
solve_ICCG2 （ 3/3 ） : METHOD=2
!C
!C***************************************************** ITERATION ITR= N
do L= 1, ITR
!C
!C +---+
!C | {z}= [Minv]{r} |
!C +---+
!C===
do i= 1, N
W(i,Z)= W(i,R) enddo
do i= 1, N WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)
WVAL= WVAL - ALlu0(k) * W(itemL(k),Z) enddo
W(i,Z)= WVAL * Dlu0(i) enddo
do i= N, 1, -1 SW = 0.0d0
do k= indexU(i-1)+1, indexU(i) SW= SW + AUlu0(k) * W(itemU(k),Z) enddo
W(i,Z)= W(i,Z) - Dlu0(i)*SW enddo
!C===
これ以下の処理は「
solve_ICCG
」と全く同じCompute r
(0)= b-[A]x
(0)for i= 1, 2, …
solve [M]z
(i-1)= r
(i-1)
i-1= r
(i-1)z
(i-1)if i=1
p
(1)= z
(0)else
i-1=
i-1/
i-2p
(i)= z
(i-1)+
i-1p
(i-1)endif
q
(i)= [A]p
(i)
i=
i-1/p
(i)q
(i)x
(i)= x
(i-1)+
ip
(i)r
(i)= r
(i-1)-
iq
(i)check convergence |r|
end
solve_ICCG2 （ 3/3 ） : METHOD=2
!C
!C************************************************* ITERATION ITR= N
do L= 1, ITR
!C
!C +---+
!C | {z}= [Minv]{r} |
!C +---+
!C===
do i= 1, N
W(i,Z)= W(i,R) enddo
do i= 1, N WVAL= W(i,Z)
do k= indexL(i-1)+1, indexL(i)
WVAL= WVAL - ALlu0(k) * W(itemL(k),Z) enddo
W(i,Z)= WVAL * Dlu0(i) enddo
do i= N, 1, -1 SW = 0.0d0
do k= indexU(i-1)+1, indexU(i) SW= SW + AUlu0(k) * W(itemU(k),Z) enddo
W(i,Z)= W(i,Z) - Dlu0(i)*SW enddo
!C===
実は，ALlu0，AUlu0の値はAL，AUと全く同じである。
METHOD=1
，METHOD=2
の答え（反復回数）は同じ M z LDL T z r
L z r
DL T z z
前進代入
Forward Substitution
後退代入
Backward Substitution
不完全修正コレスキー分解 現在のようなメッシュの場合
jk k j
k ik ij
ij
a l d l
l
1
1
n i 1 , 2 , ,
1 , , 2 ,
1
i
j
1 1 1
1
2
i k iik ik ii
i
a l d l
d
j i
j
□を満たすような（